home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / C / Applications / POV-Ray 3.0.2 / src / SOURCE / SOR.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-10-09  |  25.5 KB  |  1,312 lines  |  [TEXT/CWIE]
  1. /****************************************************************************
  2. *                   sor.c
  3. *
  4. *  This module implements functions that manipulate surfaces of revolution.
  5. *
  6. *  This module was written by Dieter Bayer [DB].
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /****************************************************************************
  27. *
  28. *  Explanation:
  29. *
  30. *    The surface of revolution primitive is defined by a set of points
  31. *    in 2d space wich are interpolated by cubic splines. The resulting
  32. *    2d function is rotated about an axis to form the final object.
  33. *
  34. *    All calculations are done in the object's (u,v,w)-coordinate system
  35. *    with the (w)-axis being the rotation axis.
  36. *
  37. *    One spline segment in the (r,w)-plane is given by the equation
  38. *
  39. *      r^2 = f(w) = A * w * w * w + B * w * w + C * w + D.
  40. *
  41. *    To intersect a ray R = P + k * D transformed into the object's
  42. *    coordinate system with the surface of revolution, the equation
  43. *
  44. *      (Pu + k * Du)^2 + (Pv + k * Dv)^2 = f(Pw + k * Dw)
  45. *
  46. *    has to be solved for k (cubic polynomial in k).
  47. *
  48. *    Note that Pu, Pv, Pw and Du, Dv, Dw denote the coordinates
  49. *    of the vectors P and D.
  50. *
  51. *  Syntax:
  52. *
  53. *    revolution
  54. *    {
  55. *      number_of_points,
  56. *
  57. *      <P[0]>, <P[1]>, ..., <P[n-1]>
  58. *
  59. *      [ open ]
  60. *    }
  61. *
  62. *    Note that the P[i] are 2d vectors where u corresponds to the radius
  63. *    and v to the height.
  64. *
  65. *    Note that the first and last point, i.e. P[0] and P[n-1], are used
  66. *    to determine the derivatives at the end point.
  67. *
  68. *    Note that the x coordinate of a point corresponds to the radius and
  69. *    the y coordinate to the height; the z coordinate isn't used.
  70. *
  71. *  ---
  72. *
  73. *  Ideas for the surface of revolution were taken from:
  74. *
  75. *    P. Burger and D. Gillies, "Rapid Ray Tracing of General Surfaces
  76. *    of Revolution", New Advances in Computer Graphics, Proceedings
  77. *    of CG International '89, R. A. Earnshaw, B. Wyvill (Eds.),
  78. *    Springer, ..., pp. 523-531
  79. *
  80. *  ---
  81. *
  82. *  May 1994 : Creation. [DB]
  83. *
  84. *****************************************************************************/
  85.  
  86. #include "frame.h"
  87. #include "povray.h"
  88. #include "vector.h"
  89. #include "povproto.h"
  90. #include "bbox.h"
  91. #include "polysolv.h"
  92. #include "matrices.h"
  93. #include "objects.h"
  94. #include "sor.h"
  95.  
  96.  
  97.  
  98. /*****************************************************************************
  99. * Local preprocessor defines
  100. ******************************************************************************/
  101.  
  102. /* Minimal intersection depth for a valid intersection. */
  103.  
  104. #define DEPTH_TOLERANCE 1.0e-4
  105.  
  106. /* |Coefficients| < COEFF_LIMIT are regarded to be 0. */
  107.  
  108. #define COEFF_LIMIT 1.0e-20
  109.  
  110. /* Part of the surface of revolution hit. */
  111.  
  112. #define BASE_PLANE 1
  113. #define CAP_PLANE  2
  114. #define CURVE      3
  115.  
  116. /* Max. number of intersecions per spline segment. */
  117.  
  118. #define MAX_INTERSECTIONS_PER_SEGMENT 4
  119.  
  120.  
  121.  
  122. /*****************************************************************************
  123. * Local typedefs
  124. ******************************************************************************/
  125.  
  126.  
  127.  
  128. /*****************************************************************************
  129. * Static functions
  130. ******************************************************************************/
  131.  
  132. static int  intersect_sor PARAMS((RAY *Ray, SOR *Sor, ISTACK *Depth_Stack));
  133. static int  All_Sor_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  134. static int  Inside_Sor PARAMS((VECTOR point, OBJECT *Object));
  135. static void Sor_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  136. static void *Copy_Sor PARAMS((OBJECT *Object));
  137. static void Translate_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  138. static void Rotate_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  139. static void Scale_Sor PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  140. static void Transform_Sor PARAMS((OBJECT *Object, TRANSFORM *Trans));
  141. static void Invert_Sor PARAMS((OBJECT *Object));
  142. static void Destroy_Sor PARAMS((OBJECT *Object));
  143.  
  144. static int test_hit PARAMS((SOR *, RAY *, ISTACK *, DBL, int, int));
  145.  
  146. /*****************************************************************************
  147. * Local variables
  148. ******************************************************************************/
  149.  
  150. static METHODS Sor_Methods =
  151. {
  152.   All_Sor_Intersections,
  153.   Inside_Sor, Sor_Normal,
  154.   Copy_Sor,
  155.   Translate_Sor, Rotate_Sor,
  156.   Scale_Sor, Transform_Sor, Invert_Sor, Destroy_Sor
  157. };
  158.  
  159.  
  160.  
  161. /*****************************************************************************
  162. *
  163. * FUNCTION
  164. *
  165. *   All_Sor_Intersections
  166. *
  167. * INPUT
  168. *
  169. *   Object      - Object
  170. *   Ray         - Ray
  171. *   Depth_Stack - Intersection stack
  172. *
  173. * OUTPUT
  174. *
  175. *   Depth_Stack
  176. *   
  177. * RETURNS
  178. *
  179. *   int - TRUE, if a intersection was found
  180. *   
  181. * AUTHOR
  182. *
  183. *   Dieter Bayer
  184. *   
  185. * DESCRIPTION
  186. *
  187. *   Determine ray/surface of revolution intersection and
  188. *   clip intersection found.
  189. *
  190. * CHANGES
  191. *
  192. *   May 1994 : Creation.
  193. *   Oct 1996 : Changed code to include faster version. [DB]
  194. *
  195. ******************************************************************************/
  196.  
  197. static int All_Sor_Intersections(Object, Ray, Depth_Stack)
  198. OBJECT *Object;
  199. RAY *Ray;
  200. ISTACK *Depth_Stack;
  201. {
  202.   Increase_Counter(stats[Ray_Sor_Tests]);
  203.  
  204.   if (intersect_sor(Ray, (SOR *)Object, Depth_Stack))
  205.   {
  206.     Increase_Counter(stats[Ray_Sor_Tests_Succeeded]);
  207.  
  208.     return(TRUE);
  209.   }
  210.  
  211.   return(FALSE);
  212. }
  213.  
  214.  
  215.  
  216. /*****************************************************************************
  217. *
  218. * FUNCTION
  219. *
  220. *   intersect_sor
  221. *
  222. * INPUT
  223. *
  224. *   Ray          - Ray
  225. *   Sor   - Sor
  226. *   Intersection - Sor intersection structure
  227. *   
  228. * OUTPUT
  229. *
  230. *   Intersection
  231. *   
  232. * RETURNS
  233. *
  234. *   int - Number of intersections found
  235. *
  236. * AUTHOR
  237. *
  238. *   Dieter Bayer
  239. *
  240. * DESCRIPTION
  241. *
  242. *   Determine ray/surface of revolution intersection.
  243. *
  244. *   NOTE: The curve is rotated about the y-axis!
  245. *         Order reduction cannot be used.
  246. *
  247. * CHANGES
  248. *
  249. *   May 1994 : Creation.
  250. *   Oct 1996 : Changed code to include faster version. [DB]
  251. *
  252. ******************************************************************************/
  253.  
  254. static int intersect_sor(Ray, Sor, Depth_Stack)
  255. RAY *Ray;
  256. SOR *Sor;
  257. ISTACK *Depth_Stack;
  258. {
  259.   int cnt;
  260.   int found, j, n;
  261.   DBL a, b, k, h, len, u, v, r0;
  262.   DBL x[4], y[3];
  263.   DBL best;
  264.   VECTOR P, D;
  265.   BCYL_INT *intervals;
  266.   SOR_SPLINE_ENTRY *Entry;
  267.  
  268.   /* Transform the ray into the surface of revolution space. */
  269.  
  270.   MInvTransPoint(P, Ray->Initial, Sor->Trans);
  271.  
  272.   MInvTransDirection(D, Ray->Direction, Sor->Trans);
  273.  
  274.   VLength(len, D);
  275.  
  276.   VInverseScaleEq(D, len);
  277.  
  278.   /* Test if ray misses object's bounds. */
  279.  
  280. #ifdef SOR_EXTRA_STATS
  281.   Increase_Counter(stats[Sor_Bound_Tests]);
  282. #endif
  283.  
  284.   if (((D[Y] >= 0.0) && (P[Y] >  Sor->Height2)) ||
  285.       ((D[Y] <= 0.0) && (P[Y] <  Sor->Height1)) ||
  286.       ((D[X] >= 0.0) && (P[X] >  Sor->Radius2)) ||
  287.       ((D[X] <= 0.0) && (P[X] < -Sor->Radius2)))
  288.   {
  289.     return(FALSE);
  290.   }
  291.  
  292.   /* Get distance of the ray from rotation axis (= y axis). */
  293.  
  294.   r0 = P[X] * D[Z] - P[Z] * D[X];
  295.  
  296.   if ((a = D[X] * D[X] + D[Z] * D[Z]) > 0.0)
  297.   {
  298.     r0 /= sqrt(a);
  299.   }
  300.  
  301.   /* Test if ray misses object's bounds. */
  302.  
  303.   if (r0 > Sor->Radius2)
  304.   {
  305.     return(FALSE);
  306.   }
  307.  
  308.   /* Intersect all cylindrical bounds. */
  309.  
  310.   if ((cnt = Intersect_BCyl(Sor->Spline->BCyl, P, D)) == 0)
  311.   {
  312.     return(FALSE);
  313.   }
  314.  
  315. #ifdef SOR_EXTRA_STATS
  316.   Increase_Counter(stats[Sor_Bound_Tests_Succeeded]);
  317. #endif
  318.  
  319.   /* Test base/cap plane. */
  320.  
  321.   found = FALSE;
  322.  
  323.   best = BOUND_HUGE;
  324.  
  325.   if (Test_Flag(Sor, CLOSED_FLAG) && (fabs(D[Y]) > EPSILON))
  326.   {
  327.     /* Test base plane. */
  328.  
  329.     if (Sor->Base_Radius_Squared > DEPTH_TOLERANCE)
  330.     {
  331.       k = (Sor->Height1 - P[Y]) / D[Y];
  332.  
  333.       u = P[X] + k * D[X];
  334.       v = P[Z] + k * D[Z];
  335.  
  336.       b = u * u + v * v;
  337.  
  338.       if (b <= Sor->Base_Radius_Squared)
  339.       {
  340.         if (test_hit(Sor, Ray, Depth_Stack, k / len, BASE_PLANE, 0))
  341.         {
  342.           found = TRUE;
  343.  
  344.           if (k < best)
  345.           {
  346.             best = k;
  347.           }
  348.         }
  349.       }
  350.     }
  351.  
  352.     /* Test cap plane. */
  353.  
  354.     if (Sor->Cap_Radius_Squared > DEPTH_TOLERANCE)
  355.     {
  356.       k = (Sor->Height1 - P[Y]) / D[Y];
  357.  
  358.       u = P[X] + k * D[X];
  359.       v = P[Z] + k * D[Z];
  360.  
  361.       b = u * u + v * v;
  362.  
  363.       if (b <= Sor->Cap_Radius_Squared)
  364.       {
  365.         if (test_hit(Sor, Ray, Depth_Stack, k / len, CAP_PLANE, 0))
  366.         {
  367.           found = TRUE;
  368.  
  369.           if (k < best)
  370.           {
  371.             best = k;
  372.           }
  373.         }
  374.       }
  375.     }
  376.   }
  377.  
  378.   /* Step through the list of intersections. */
  379.  
  380.   intervals = Sor->Spline->BCyl->intervals;
  381.  
  382.   for (j = 0; j < cnt; j++)
  383.   {
  384.     /* Get current segment. */
  385.  
  386.     Entry = &Sor->Spline->Entry[intervals[j].n];
  387.  
  388.     /* If we already have the best intersection we may exit. */
  389.  
  390.     if (!(Sor->Type & IS_CHILD_OBJECT) && (intervals[j].d[0] > best))
  391.     {
  392.       break;
  393.     }
  394.  
  395.     /* Cubic curve. */
  396.  
  397.     x[0] = Entry->A * D[Y] * D[Y] * D[Y];
  398.  
  399. /*
  400.     x[1] = D[Y] * D[Y] * (3.0 * Entry->A * P[Y] + Entry->B) - D[X] * D[X] - D[Z] * D[Z];
  401. */
  402.     x[1] = D[Y] * D[Y] * (3.0 * Entry->A * P[Y] + Entry->B) - a;
  403.  
  404.     x[2] = D[Y] * (P[Y] * (3.0 * Entry->A * P[Y] + 2.0 * Entry->B) + Entry->C) - 2.0 * (P[X] * D[X] + P[Z] * D[Z]);
  405.  
  406.     x[3] = P[Y] * (P[Y] * (Entry->A * P[Y] + Entry->B) + Entry->C) + Entry->D - P[X] * P[X] - P[Z] * P[Z];
  407.  
  408.     n = Solve_Polynomial(3, x, y, Test_Flag(Sor, STURM_FLAG), 0.0);
  409.  
  410.     while (n--)
  411.     {
  412.       k = y[n];
  413.  
  414.       h = P[Y] + k * D[Y];
  415.  
  416.       if ((h >= Sor->Spline->BCyl->height[Sor->Spline->BCyl->entry[intervals[j].n].h1]) &&
  417.           (h <= Sor->Spline->BCyl->height[Sor->Spline->BCyl->entry[intervals[j].n].h2]))
  418.       {
  419.         if (test_hit(Sor, Ray, Depth_Stack, k / len, CURVE, intervals[j].n))
  420.         {
  421.           found = TRUE;
  422.  
  423.           if (y[n] < best)
  424.           {
  425.             best = k;
  426.           }
  427.         }
  428.       }
  429.     }
  430.   }
  431.  
  432.   return(found);
  433. }
  434.  
  435.  
  436.  
  437. /*****************************************************************************
  438. *
  439. * FUNCTION
  440. *
  441. *   Inside_Sor
  442. *
  443. * INPUT
  444. *
  445. *   IPoint - Intersection point
  446. *   Object - Object
  447. *   
  448. * OUTPUT
  449. *   
  450. * RETURNS
  451. *
  452. *   int - TRUE if inside
  453. *   
  454. * AUTHOR
  455. *
  456. *   Dieter Bayer
  457. *   
  458. * DESCRIPTION
  459. *
  460. *   Return true if point lies inside the surface of revolution.
  461. *
  462. * CHANGES
  463. *
  464. *   May 1994 : Creation.
  465. *
  466. ******************************************************************************/
  467.  
  468. static int Inside_Sor(IPoint, Object)
  469. VECTOR IPoint;
  470. OBJECT *Object;
  471. {
  472.   int i;
  473.   DBL r0, r;
  474.   VECTOR P;
  475.   SOR *Sor = (SOR *)Object;
  476.   SOR_SPLINE_ENTRY *Entry;
  477.  
  478.   /* Transform the point into the surface of revolution space. */
  479.  
  480.   MInvTransPoint(P, IPoint, Sor->Trans);
  481.  
  482.   /* Test if we are inside the cylindrical bound. */
  483.  
  484.   if ((P[Y] >= Sor->Height1) && (P[Y] <= Sor->Height2))
  485.   {
  486.     r0 = P[X] * P[X] + P[Z] * P[Z];
  487.  
  488.     /* Test if we are inside the cylindrical bound. */
  489.  
  490.     if (r0 <= Sqr(Sor->Radius2))
  491.     {
  492.       /* Now find the segment the point is in. */
  493.  
  494.       for (i = 0; i < Sor->Number; i++)
  495.       {
  496.         Entry = &Sor->Spline->Entry[i];
  497.  
  498.         if ((P[Y] >= Sor->Spline->BCyl->height[Sor->Spline->BCyl->entry[i].h1]) &&
  499.             (P[Y] <= Sor->Spline->BCyl->height[Sor->Spline->BCyl->entry[i].h2]))
  500.         {
  501.           break;
  502.         }
  503.       }
  504.  
  505.       /* Have we found any segment? */
  506.  
  507.       if (i < Sor->Number)
  508.       {
  509.         r = P[Y] * (P[Y] * (P[Y] * Entry->A + Entry->B) + Entry->C) + Entry->D;
  510.  
  511.         if (r0 <= r)
  512.         {
  513.           /* We're inside. */
  514.  
  515.           return(!Test_Flag(Sor, INVERTED_FLAG));
  516.         }
  517.       }
  518.     }
  519.   }
  520.  
  521.   /* We're outside. */
  522.  
  523.   return(Test_Flag(Sor, INVERTED_FLAG));
  524. }
  525.  
  526.  
  527.  
  528. /*****************************************************************************
  529. *
  530. * FUNCTION
  531. *
  532. *   Sor_Normal
  533. *
  534. * INPUT
  535. *
  536. *   Result - Normal vector
  537. *   Object - Object
  538. *   Inter  - Intersection found
  539. *   
  540. * OUTPUT
  541. *
  542. *   Result
  543. *
  544. * RETURNS
  545. *   
  546. * AUTHOR
  547. *
  548. *   Dieter Bayer
  549. *   
  550. * DESCRIPTION
  551. *
  552. *   Calculate the normal of the surface of revolution in a given point.
  553. *
  554. * CHANGES
  555. *
  556. *   May 1994 : Creation.
  557. *
  558. ******************************************************************************/
  559.  
  560. static void Sor_Normal(Result, Object, Inter)
  561. OBJECT *Object;
  562. VECTOR Result;
  563. INTERSECTION *Inter;
  564. {
  565.   DBL k;
  566.   VECTOR P;
  567.   SOR *Sor = (SOR *)Object;
  568.   SOR_SPLINE_ENTRY *Entry;
  569.   VECTOR N;
  570.  
  571.   switch (Inter->i1)
  572.   {
  573.     case CURVE:
  574.  
  575.       /* Transform the intersection point into the surface of revolution space. */
  576.  
  577.       MInvTransPoint(P, Inter->IPoint, Sor->Trans);
  578.  
  579.       if (P[X] * P[X] + P[Z] * P[Z] > DEPTH_TOLERANCE)
  580.       {
  581.         Entry = &Sor->Spline->Entry[Inter->i2];
  582.  
  583.         k = 0.5 * (P[Y] * (3.0 * Entry->A * P[Y] + 2.0 * Entry->B) + Entry->C);
  584.  
  585.         N[X] = P[X];
  586.         N[Y] = -k;
  587.         N[Z] = P[Z];
  588.       }
  589.  
  590.       break;
  591.  
  592.     case BASE_PLANE:
  593.  
  594.       Make_Vector(N, 0.0, -1.0, 0.0);
  595.  
  596.       break;
  597.  
  598.  
  599.     case CAP_PLANE:
  600.  
  601.       Make_Vector(N, 0.0, 1.0, 0.0);
  602.  
  603.       break;
  604.   }
  605.  
  606.   /* Transform the normal out of the surface of revolution space. */
  607.  
  608.   MTransNormal(Result, N, Sor->Trans);
  609.  
  610.   VNormalize(Result, Result);
  611. }
  612.  
  613.  
  614.  
  615. /*****************************************************************************
  616. *
  617. * FUNCTION
  618. *
  619. *   Translate_Sor
  620. *
  621. * INPUT
  622. *
  623. *   Object - Object
  624. *   Vector - Translation vector
  625. *   
  626. * OUTPUT
  627. *
  628. *   Object
  629. *   
  630. * RETURNS
  631. *   
  632. * AUTHOR
  633. *
  634. *   Dieter Bayer
  635. *   
  636. * DESCRIPTION
  637. *
  638. *   Translate a surface of revolution.
  639. *
  640. * CHANGES
  641. *
  642. *   May 1994 : Creation.
  643. *
  644. ******************************************************************************/
  645.  
  646. static void Translate_Sor(Object, Vector, Trans)
  647. OBJECT *Object;
  648. VECTOR Vector;
  649. TRANSFORM *Trans;
  650. {
  651.   Transform_Sor(Object, Trans);
  652. }
  653.  
  654.  
  655.  
  656. /*****************************************************************************
  657. *
  658. * FUNCTION
  659. *
  660. *   Rotate_Sor
  661. *
  662. * INPUT
  663. *
  664. *   Object - Object
  665. *   Vector - Rotation vector
  666. *   
  667. * OUTPUT
  668. *
  669. *   Object
  670. *   
  671. * RETURNS
  672. *   
  673. * AUTHOR
  674. *
  675. *   Dieter Bayer
  676. *   
  677. * DESCRIPTION
  678. *
  679. *   Rotate a surface of revolution.
  680. *
  681. * CHANGES
  682. *
  683. *   May 1994 : Creation.
  684. *
  685. ******************************************************************************/
  686.  
  687. static void Rotate_Sor(Object, Vector, Trans)
  688. OBJECT *Object;
  689. VECTOR Vector;
  690. TRANSFORM *Trans;
  691. {
  692.   Transform_Sor(Object, Trans);
  693. }
  694.  
  695.  
  696.  
  697. /*****************************************************************************
  698. *
  699. * FUNCTION
  700. *
  701. *   Scale_Sor
  702. *
  703. * INPUT
  704. *
  705. *   Object - Object
  706. *   Vector - Scaling vector
  707. *   
  708. * OUTPUT
  709. *
  710. *   Object
  711. *   
  712. * RETURNS
  713. *   
  714. * AUTHOR
  715. *
  716. *   Dieter Bayer
  717. *   
  718. * DESCRIPTION
  719. *
  720. *   Scale a surface of revolution.
  721. *
  722. * CHANGES
  723. *
  724. *   May 1994 : Creation.
  725. *
  726. ******************************************************************************/
  727.  
  728. static void Scale_Sor(Object, Vector, Trans)
  729. OBJECT *Object;
  730. VECTOR Vector;
  731. TRANSFORM *Trans;
  732. {
  733.   Transform_Sor(Object, Trans);
  734. }
  735.  
  736.  
  737.  
  738. /*****************************************************************************
  739. *
  740. * FUNCTION
  741. *
  742. *   Transform_Sor
  743. *
  744. * INPUT
  745. *
  746. *   Object - Object
  747. *   Trans  - Transformation to apply
  748. *   
  749. * OUTPUT
  750. *
  751. *   Object
  752. *   
  753. * RETURNS
  754. *   
  755. * AUTHOR
  756. *
  757. *   Dieter Bayer
  758. *   
  759. * DESCRIPTION
  760. *
  761. *   Transform a surface of revolution and recalculate its bounding box.
  762. *
  763. * CHANGES
  764. *
  765. *   May 1994 : Creation.
  766. *
  767. ******************************************************************************/
  768.  
  769. static void Transform_Sor(Object, Trans)
  770. OBJECT *Object;
  771. TRANSFORM *Trans;
  772. {
  773.   Compose_Transforms(((SOR *)Object)->Trans, Trans);
  774.  
  775.   Compute_Sor_BBox((SOR *)Object);
  776. }
  777.  
  778.  
  779.  
  780. /*****************************************************************************
  781. *
  782. * FUNCTION
  783. *
  784. *   Invert_Sor
  785. *
  786. * INPUT
  787. *
  788. *   Object - Object
  789. *   
  790. * OUTPUT
  791. *
  792. *   Object
  793. *   
  794. * RETURNS
  795. *   
  796. * AUTHOR
  797. *
  798. *   Dieter Bayer
  799. *   
  800. * DESCRIPTION
  801. *
  802. *   Invert a surface of revolution.
  803. *
  804. * CHANGES
  805. *
  806. *   May 1994 : Creation.
  807. *
  808. ******************************************************************************/
  809.  
  810. static void Invert_Sor(Object)
  811. OBJECT *Object;
  812. {
  813.   Invert_Flag(Object, INVERTED_FLAG);
  814. }
  815.  
  816.  
  817.  
  818. /*****************************************************************************
  819. *
  820. * FUNCTION
  821. *
  822. *   Create_Sor
  823. *
  824. * INPUT
  825. *   
  826. * OUTPUT
  827. *   
  828. * RETURNS
  829. *
  830. *   SOR * - new surface of revolution
  831. *   
  832. * AUTHOR
  833. *
  834. *   Dieter Bayer
  835. *   
  836. * DESCRIPTION
  837. *
  838. *   Create a new surface of revolution.
  839. *
  840. * CHANGES
  841. *
  842. *   May 1994 : Creation.
  843. *
  844. ******************************************************************************/
  845.  
  846. SOR *Create_Sor()
  847. {
  848.   SOR *New;
  849.  
  850.   New = (SOR *)POV_MALLOC(sizeof(SOR), "surface of revolution");
  851.  
  852.   INIT_OBJECT_FIELDS(New,SOR_OBJECT,&Sor_Methods)
  853.  
  854.   New->Trans = Create_Transform();
  855.  
  856.   New->Spline = NULL;
  857.  
  858.   New->Radius2             =
  859.   New->Base_Radius_Squared =
  860.   New->Cap_Radius_Squared  = 0.0;
  861.  
  862.   return(New);
  863. }
  864.  
  865.  
  866.  
  867. /*****************************************************************************
  868. *
  869. * FUNCTION
  870. *
  871. *   Copy_Sor
  872. *
  873. * INPUT
  874. *
  875. *   Object - Object
  876. *   
  877. * OUTPUT
  878. *   
  879. * RETURNS
  880. *
  881. *   void * - New surface of revolution
  882. *   
  883. * AUTHOR
  884. *
  885. *   Dieter Bayer
  886. *   
  887. * DESCRIPTION
  888. *
  889. *   Copy a surface of revolution structure.
  890. *
  891. *   NOTE: The splines are not copied, only the number of references is
  892. *         counted, so that Destray_Sor() knows if they can be destroyed.
  893. *
  894. * CHANGES
  895. *
  896. *   May 1994 : Creation.
  897. *
  898. *   Sep 1994 : fixed memory leakage [DB]
  899. *
  900. ******************************************************************************/
  901.  
  902. static void *Copy_Sor(Object)
  903. OBJECT *Object;
  904. {
  905.   SOR *New, *Sor = (SOR *)Object;
  906.  
  907.   New = Create_Sor();
  908.  
  909.   /* Get rid of the transformation created in Create_Sor(). */
  910.  
  911.   Destroy_Transform(New->Trans);
  912.  
  913.   /* Copy surface of revolution. */
  914.  
  915.   *New = *Sor;
  916.  
  917.   New->Trans = Copy_Transform(Sor->Trans);
  918.  
  919.   New->Spline->References++;
  920.  
  921.   return(New);
  922. }
  923.  
  924.  
  925.  
  926. /*****************************************************************************
  927. *
  928. * FUNCTION
  929. *
  930. *   Destroy_Sor
  931. *
  932. * INPUT
  933. *
  934. *   Object - Object
  935. *   
  936. * OUTPUT
  937. *
  938. *   Object
  939. *   
  940. * RETURNS
  941. *   
  942. * AUTHOR
  943. *
  944. *   Dieter Bayer
  945. *   
  946. * DESCRIPTION
  947. *
  948. *   Destroy a surface of revolution.
  949. *
  950. *   NOTE: The splines are destroyed if they are no longer used by any copy.
  951. *
  952. * CHANGES
  953. *
  954. *   May 1994 : Creation.
  955. *
  956. ******************************************************************************/
  957.  
  958. static void Destroy_Sor (Object)
  959. OBJECT *Object;
  960. {
  961.   SOR *Sor = (SOR *)Object;
  962.  
  963.   Destroy_Transform(Sor->Trans);
  964.  
  965.   if (--(Sor->Spline->References) == 0)
  966.   {
  967.     Destroy_BCyl(Sor->Spline->BCyl);
  968.  
  969.     POV_FREE(Sor->Spline->Entry);
  970.  
  971.     POV_FREE(Sor->Spline);
  972.   }
  973.  
  974.   POV_FREE(Object);
  975. }
  976.  
  977.  
  978.  
  979. /*****************************************************************************
  980. *
  981. * FUNCTION
  982. *
  983. *   Compute_Sor_BBox
  984. *
  985. * INPUT
  986. *
  987. *   Sor - Sor
  988. *   
  989. * OUTPUT
  990. *
  991. *   Sor
  992. *   
  993. * RETURNS
  994. *   
  995. * AUTHOR
  996. *
  997. *   Dieter Bayer
  998. *   
  999. * DESCRIPTION
  1000. *
  1001. *   Calculate the bounding box of a surface of revolution.
  1002. *
  1003. * CHANGES
  1004. *
  1005. *   May 1994 : Creation.
  1006. *
  1007. ******************************************************************************/
  1008.  
  1009. void Compute_Sor_BBox(Sor)
  1010. SOR *Sor;
  1011. {
  1012.   Make_BBox(Sor->BBox, -Sor->Radius2, Sor->Height1, -Sor->Radius2,
  1013.     2.0 * Sor->Radius2, Sor->Height2 - Sor->Height1, 2.0 * Sor->Radius2);
  1014.  
  1015.   Recompute_BBox(&Sor->BBox, Sor->Trans);
  1016. }
  1017.  
  1018.  
  1019.  
  1020. /*****************************************************************************
  1021. *
  1022. * FUNCTION
  1023. *
  1024. *   Compute_Sor
  1025. *
  1026. * INPUT
  1027. *
  1028. *   Sor - Sor
  1029. *   P          - Points defining surface of revolution
  1030. *   
  1031. * OUTPUT
  1032. *
  1033. *   Sor
  1034. *   
  1035. * RETURNS
  1036. *   
  1037. * AUTHOR
  1038. *
  1039. *   Dieter Bayer, June 1994
  1040. *   
  1041. * DESCRIPTION
  1042. *
  1043. *   Calculate the spline segments of a surface of revolution
  1044. *   from a set of points.
  1045. *
  1046. *   Note that the number of points in the surface of revolution has to be set.
  1047. *
  1048. * CHANGES
  1049. *
  1050. *   May 1994 : Creation.
  1051. *
  1052. ******************************************************************************/
  1053.  
  1054. void Compute_Sor(Sor, P)
  1055. SOR *Sor;
  1056. UV_VECT *P;
  1057. {
  1058.   int i, n;
  1059.   DBL *tmp_r1;
  1060.   DBL *tmp_r2;
  1061.   DBL *tmp_h1;
  1062.   DBL *tmp_h2;
  1063.   DBL A, B, C, D, w, k[4];
  1064.   DBL x[4], xmax, xmin;
  1065.   DBL y[2], ymax, ymin;
  1066.   DBL c[3], r[2];
  1067.   MATRIX Mat;
  1068.  
  1069.   /* Allocate Sor->Number segments. */
  1070.  
  1071.   if (Sor->Spline == NULL)
  1072.   {
  1073.     Sor->Spline = (SOR_SPLINE *)POV_MALLOC(sizeof(SOR_SPLINE), "spline segments of surface of revoluion");
  1074.  
  1075.     Sor->Spline->References = 1;
  1076.  
  1077.     Sor->Spline->Entry = (SOR_SPLINE_ENTRY *)POV_MALLOC(Sor->Number*sizeof(SOR_SPLINE_ENTRY), "spline segments of surface of revoluion");
  1078.   }
  1079.   else
  1080.   {
  1081.     Error("Surface of revolution segments are already defined.\n");
  1082.   }
  1083.  
  1084.   /* Allocate temporary lists. */
  1085.  
  1086.   tmp_r1 = (DBL *)POV_MALLOC(Sor->Number * sizeof(DBL), "temp lathe data");
  1087.   tmp_r2 = (DBL *)POV_MALLOC(Sor->Number * sizeof(DBL), "temp lathe data");
  1088.   tmp_h1 = (DBL *)POV_MALLOC(Sor->Number * sizeof(DBL), "temp lathe data");
  1089.   tmp_h2 = (DBL *)POV_MALLOC(Sor->Number * sizeof(DBL), "temp lathe data");
  1090.  
  1091.   /* We want to know the size of the overall bounding cylinder. */
  1092.  
  1093.   xmax = ymax = -BOUND_HUGE;
  1094.   xmin = ymin =  BOUND_HUGE;
  1095.  
  1096.   /* Calculate segments, i.e. cubic patches. */
  1097.  
  1098.   for (i = 0; i < Sor->Number; i++)
  1099.   {
  1100.     if ((fabs(P[i+2][Y] - P[i][Y]) < EPSILON) ||
  1101.         (fabs(P[i+3][Y] - P[i+1][Y]) < EPSILON))
  1102.     {
  1103.       Error("Incorrect point in surface of revolution.\n");
  1104.     }
  1105.  
  1106.     /* Use cubic interpolation. */
  1107.  
  1108.     k[0] = P[i+1][X] * P[i+1][X];
  1109.     k[1] = P[i+2][X] * P[i+2][X];
  1110.     k[2] = (P[i+2][X] - P[i][X]) / (P[i+2][Y] - P[i][Y]);
  1111.     k[3] = (P[i+3][X] - P[i+1][X]) / (P[i+3][Y] - P[i+1][Y]);
  1112.  
  1113.     k[2] *= 2.0 * P[i+1][X];
  1114.     k[3] *= 2.0 * P[i+2][X];
  1115.  
  1116.     w = P[i+1][Y];
  1117.  
  1118.     Mat[0][0] = w * w * w;
  1119.     Mat[0][1] = w * w;
  1120.     Mat[0][2] = w;
  1121.     Mat[0][3] = 1.0;
  1122.  
  1123.     Mat[2][0] = 3.0 * w * w;
  1124.     Mat[2][1] = 2.0 * w;
  1125.     Mat[2][2] = 1.0;
  1126.     Mat[2][3] = 0.0;
  1127.  
  1128.     w = P[i+2][Y];
  1129.  
  1130.     Mat[1][0] = w * w * w;
  1131.     Mat[1][1] = w * w;
  1132.     Mat[1][2] = w;
  1133.     Mat[1][3] = 1.0;
  1134.  
  1135.     Mat[3][0] = 3.0 * w * w;
  1136.     Mat[3][1] = 2.0 * w;
  1137.     Mat[3][2] = 1.0;
  1138.     Mat[3][3] = 0.0;
  1139.  
  1140.     MInvers(Mat, Mat);
  1141.  
  1142.     /* Calculate coefficients of cubic patch. */
  1143.  
  1144.     A = k[0] * Mat[0][0] + k[1] * Mat[0][1] + k[2] * Mat[0][2] + k[3] * Mat[0][3];
  1145.     B = k[0] * Mat[1][0] + k[1] * Mat[1][1] + k[2] * Mat[1][2] + k[3] * Mat[1][3];
  1146.     C = k[0] * Mat[2][0] + k[1] * Mat[2][1] + k[2] * Mat[2][2] + k[3] * Mat[2][3];
  1147.     D = k[0] * Mat[3][0] + k[1] * Mat[3][1] + k[2] * Mat[3][2] + k[3] * Mat[3][3];
  1148.  
  1149.     if (fabs(A) < EPSILON) A = 0.0;
  1150.     if (fabs(B) < EPSILON) B = 0.0;
  1151.     if (fabs(C) < EPSILON) C = 0.0;
  1152.     if (fabs(D) < EPSILON) D = 0.0;
  1153.  
  1154.     Sor->Spline->Entry[i].A = A;
  1155.     Sor->Spline->Entry[i].B = B;
  1156.     Sor->Spline->Entry[i].C = C;
  1157.     Sor->Spline->Entry[i].D = D;
  1158.  
  1159.     /* Get minimum and maximum radius**2 in current segment. */
  1160.  
  1161.     y[0] = P[i+1][Y];
  1162.     y[1] = P[i+2][Y];
  1163.  
  1164.     x[0] = x[2] = P[i+1][X];
  1165.     x[1] = x[3] = P[i+2][X];
  1166.  
  1167.     c[0] = 3.0 * A;
  1168.     c[1] = 2.0 * B;
  1169.     c[2] = C;
  1170.  
  1171.     n = Solve_Polynomial(2, c, r, FALSE, 0.0);
  1172.  
  1173.     while (n--)
  1174.     {
  1175.       if ((r[n] >= y[0]) && (r[n] <= y[1]))
  1176.       {
  1177.         x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);
  1178.       }
  1179.     }
  1180.  
  1181.     /* Set current segment's bounding cylinder. */
  1182.  
  1183.     tmp_r1[i] = min(min(x[0], x[1]), min(x[2], x[3]));
  1184.     tmp_r2[i] = max(max(x[0], x[1]), max(x[2], x[3]));
  1185.  
  1186.     tmp_h1[i] = y[0];
  1187.     tmp_h2[i] = y[1];
  1188.  
  1189.     /* Keep track of overall bounding cylinder. */
  1190.  
  1191.     xmin = min(xmin, tmp_r1[i]);
  1192.     xmax = max(xmax, tmp_r2[i]);
  1193.  
  1194.     ymin = min(ymin, tmp_h1[i]);
  1195.     ymax = max(ymax, tmp_h2[i]);
  1196.  
  1197. /*
  1198.     fprintf(stderr, "bound spline segment %d: ", i);
  1199.     fprintf(stderr, "r = %f - %f, h = %f - %f\n", tmp_r1[i], tmp_r2[i], tmp_h1[i], tmp_h2[i]);
  1200. */
  1201.   }
  1202.  
  1203.   /* Set overall bounding cylinder. */
  1204.  
  1205.   Sor->Radius1 = ymax;
  1206.   Sor->Radius2 = xmax;
  1207.  
  1208.   Sor->Height1 = ymin;
  1209.   Sor->Height2 = ymax;
  1210.  
  1211.   /* Get cap radius. */
  1212.  
  1213.   w = tmp_h1[Sor->Number-1];
  1214.  
  1215.   A = Sor->Spline->Entry[Sor->Number-1].A;
  1216.   B = Sor->Spline->Entry[Sor->Number-1].B;
  1217.   C = Sor->Spline->Entry[Sor->Number-1].C;
  1218.   D = Sor->Spline->Entry[Sor->Number-1].D;
  1219.  
  1220.   if ((Sor->Cap_Radius_Squared = w * (w * (A * w + B) + C) + D) < 0.0)
  1221.   {
  1222.     Sor->Cap_Radius_Squared = 0.0;
  1223.   }
  1224.  
  1225.   /* Get base radius. */
  1226.  
  1227.   w = tmp_h1[0];
  1228.  
  1229.   A = Sor->Spline->Entry[0].A;
  1230.   B = Sor->Spline->Entry[0].B;
  1231.   C = Sor->Spline->Entry[0].C;
  1232.   D = Sor->Spline->Entry[0].D;
  1233.  
  1234.   if ((Sor->Base_Radius_Squared = w * (w * (A * w + B) + C) + D) < 0.0)
  1235.   {
  1236.     Sor->Base_Radius_Squared = 0.0;
  1237.   }
  1238.  
  1239.   /* Get bounding cylinder. */
  1240.  
  1241.   Sor->Spline->BCyl = Create_BCyl(Sor->Number, tmp_r1, tmp_r2, tmp_h1, tmp_h2);
  1242.  
  1243.   /* Get rid of temp. memory. */
  1244.  
  1245.   POV_FREE(tmp_h2);
  1246.   POV_FREE(tmp_h1);
  1247.   POV_FREE(tmp_r2);
  1248.   POV_FREE(tmp_r1);
  1249. }
  1250.  
  1251.  
  1252.  
  1253. /*****************************************************************************
  1254. *
  1255. * FUNCTION
  1256. *
  1257. *   test_hit
  1258. *
  1259. * INPUT
  1260. *
  1261. *   Sor         - Pointer to lathe structure
  1262. *   Ray         - Current ray
  1263. *   Depth_Stack - Current depth stack
  1264. *   d, t, n     - Intersection depth, type and number
  1265. *
  1266. * OUTPUT
  1267. *
  1268. *   Depth_Stack
  1269. *
  1270. * RETURNS
  1271. *
  1272. * AUTHOR
  1273. *
  1274. *   Dieter Bayer
  1275. *
  1276. * DESCRIPTION
  1277. *
  1278. *   Test if a hit is valid and push if on the intersection depth.
  1279. *
  1280. * CHANGES
  1281. *
  1282. *   Oct 1996 : Creation.
  1283. *
  1284. ******************************************************************************/
  1285.  
  1286. static int test_hit(Sor, Ray, Depth_Stack, d, t, n)
  1287. SOR *Sor;
  1288. RAY *Ray;
  1289. ISTACK *Depth_Stack;
  1290. DBL d;
  1291. int t, n;
  1292. {
  1293.   VECTOR IPoint;
  1294.  
  1295.   if ((d > DEPTH_TOLERANCE) && (d < Max_Distance))
  1296.   {
  1297.     VEvaluateRay(IPoint, Ray->Initial, d, Ray->Direction);
  1298.  
  1299.     if (Point_In_Clip(IPoint, ((OBJECT *)Sor)->Clip))
  1300.     {
  1301.       push_entry_i1_i2(d, IPoint, (OBJECT *)Sor, t, n, Depth_Stack);
  1302.  
  1303.       return(TRUE);
  1304.     }
  1305.   }
  1306.  
  1307.   return(FALSE);
  1308. }
  1309.  
  1310.  
  1311.  
  1312.